fdr和pvalue显著均为0.05
setwd("E:\\5hmc_file\\2_5hmc_yjp_bam\\ASM\\bayes_pvalue_beta0")
group1=c("X2B_X1T","M8_M7","M6_M5","M2_M1","M48_M47","M28_M27","M30_M29","M26_M25","M35_M36","M18_M17","M20_M19","M22_M21","M40_M39","M50_M49")
sel1=paste(group1,".bayes_p.txt",sep="")

for(i in 2:length(sel1)){
  file=read.table(sel1[i],head=T,sep="\t")
  file$unitID=paste(file$chrom,file$position,sep=":")
  file$normal_bayes_fdr=p.adjust(file$normal_bayes_pvalue,method = "BH")
  file$tumor_bayes_fdr=p.adjust(file$tumor_bayes_pvalue,method = "BH")
  normal_sig_tmp1=file[file$normal_bayes_pvalue<0.05,]
  tumor_sig_tmp1=file[file$tumor_bayes_pvalue<0.05,]
  normal_sig_tmp1=data.frame(unitID=normal_sig_tmp1$unitID)
  normal_sig_tmp1$source=unlist(strsplit(gsub(".bayes_p.txt","",sel1[i]),"_"))[1]
  tumor_sig_tmp1=data.frame(unitID=tumor_sig_tmp1$unitID)
  tumor_sig_tmp1$source=unlist(strsplit(gsub(".bayes_p.txt","",sel1[i]),"_"))[2]
  tmp1=merge(normal_sig_tmp1,tumor_sig_tmp1,by="unitID",all = T)
  tmp11=merge(tmp11,tmp1,by="unitID",all=T)
  
  normal_sig_tmp2=file[file$normal_bayes_fdr<0.05,]
  tumor_sig_tmp2=file[file$tumor_bayes_fdr<0.05,]
  normal_sig_tmp2=data.frame(unitID=normal_sig_tmp2$unitID)
  normal_sig_tmp2$source=unlist(strsplit(gsub(".bayes_p.txt","",sel1[i]),"_"))[1]
  tumor_sig_tmp2=data.frame(unitID=tumor_sig_tmp2$unitID)
  tumor_sig_tmp2$source=unlist(strsplit(gsub(".bayes_p.txt","",sel1[i]),"_"))[2]
  tmp2=merge(normal_sig_tmp2,tumor_sig_tmp2,by="unitID",all = T)
  tmp22=merge(tmp22,tmp2,by="unitID",all=T)
}
header=c("unitID",paste("source",1:28,sep=""))
names(tmp11)=header
names(tmp22)=header
write.table(tmp11,"pvalue_sig.txt",quote=F,row.names = F,sep="\t")
write.table(tmp22,"fdr_sig.txt",quote=F,row.names = F,sep="\t")




setwd("E:\\5hmc_file\\2_5hmc_yjp_bam\\ASM\\bayes_pvalue_beta0")
library(tidyr)
group1=c("X2B_X1T","M8_M7","M6_M5","M2_M1","M48_M47","M28_M27","M30_M29","M26_M25","M35_M36","M18_M17","M20_M19","M22_M21","M40_M39","M50_M49")

file=read.table("pvalue_sig.txt",head=T,sep="\t")
file$num1=rowSums(!is.na(file[,2:29]),na.rm = T)
file1=unite(file,tag,source1:source28,sep="_",na.rm = T)
file1$twins=0
file1$single=0
sel=which(file1$num1>1)
file1[file1$num1==1,]$single=1


for(i in sel){
x=0
j=1
while (j <length(unlist(strsplit(file1$tag[i],"_")))) {
  str1=paste(unlist(strsplit(file1$tag[i],"_"))[j],unlist(strsplit(file1$tag[i],"_"))[j+1],sep="_")
  a=length(intersect(group1,str1))
  if(abs(a)==1){
    x=x+1
    j=j+2
  }else{
    j=j+1
    }
}
file1[i,4]=x
file1[i,5]=file1[i,]$num1-2*x
}
file1$tag2=paste(file1$twins,"twins","_",file1$single,"single",sep="")
test=data.frame(table(file1$tag2))
write.csv(test,"pvalue_sig_in_twinsORsingle.csv",quote = F,row.names = F)

###旧的循环是通过相邻两个字符串中的数字相减是否等于1来判断是否为同一对twins
for(i in sel){
x=0
j=1
while (j <length(unlist(strsplit(file1$tag[i],"_")))) {
  a=as.numeric(parse_number(unlist(strsplit(file1$tag[i],"_"))[j]))-as.numeric(parse_number(unlist(strsplit(file1$tag[i],"_"))[j+1]))
  if(abs(a)==1){
    x=x+1
    j=j+2
  }else{
    j=j+1
    }
}
file1[i,4]=x
file1[i,5]=file1[i,]$num1-2*x
}